Boundary singularities and boundary conditions for the Fokker-Planck equations 
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The boundary conditions for the Fokker-Planck equations, forward and backward ones are directly 
derived from the Chapman-Kolmogorov equation for Af -dimensional region with boundaries. The 
boundaries are assumed, in addition, to be able to absorb wandering particles or to give rise to 
fast surface transport. It is demonstrated that the boundaries break down the symmetry of random 
walks in their vicinity, leading to the boundary singularities in the corresponding kinetic coefficients. 
Eliminating these singularities we get the desired boundary conditions. As it must be the boundary 
condition for the forward Fokker-Planck equation matches the mass conservation. 
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I. THE CHAPMAN-KOLMOGOROV AND 
FOKKER-PLANCK EQUATIONS. EFFECT OF 
THE MEDIUM BOUNDARIES 

As is well-known [l], 0] Markovian stochastic processes 
are completely determined by their transition probabili- 
ties which obey the Chapman-Kolmogorov equation. The 
Kramers-Moyal expansion can be used to determine the 
Fokker-Planck equation by specifying drift vector and dif- 
fusion tensor based on the assumption of vanishing higher 
order Kramers-Moyal coefficients. 

Usually, the Fokker-Planck equations are derived im- 
plicitly assuming that the phase space of the stochastic 
variables under consideration extends to infinity so that 
so-called natural boundary conditions apply. If stochas- 
tic processes in a finite region of phase space are con- 
sidered, boundary conditions are introduced a posteri- 
ori based on apparent physical arguments leading to the 
notion of a reflecting barrier, characterized by a vanish- 
ing normal component of the probability current, an ab- 
sorbing barrier, where the probability distribution has 
to vanish, and boundary conditions at a discontinuity, 
where probability distributions and the normal compo- 
nents of the probability current have to be continuous. 
No attempts, so far, have been made to derive the Fokker- 
Planck equation simultaneously with appropriate bound- 
ary conditions from the Chapman-Kolmogorov equation. 

It is quite evident that boundaries can strongly influ- 
ence the stochastic motion of a particle in various ways 
depending on the microscopic interactions. As an exam- 
ple we mention a boundary formed by a fast diffusion 
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layer. In such a thin layer particles are able to diffuse in 
the directions tangential to the boundary on a fast time 
scale, whereas in the bulk the particles behavior should 
accurately be described by the Fokker-Planck equation. 
The theoretical treatment of the particle diffusion re- 
quires a formulation of consistent boundary conditions 
which match the internal Fokker-Planck behavior to the 
stochastic properties of the boundary layer. 

So it could be desirable to have a technique of de- 
riving the boundary conditions applying directly to the 
manner of the region boundaries affecting stochastic pro- 
cesses. In this respect we note paper [3| devoted to the 
general description of random processes near boundaries 
causing deterministic jumps, paper 0] deriving bound- 
ary conditions for the Fokker-Planck equation describ- 
ing coupled transport of photons and electrons, a seri- 
ous of papers [1, IE 0] dealing with boundary conditions 
for the advection-diffusion problem combining the Boltz- 
mann and Fokker-Planck equations and their numerical 
implementation, and also work |§| developing diffusion 
models for molecular transport across membranes via ion 
channels and wider pores in terms of random walks af- 
fected by boundaries with complex properties. In addi- 
tion paper Q actually constructs the absorbing boundary 
as a limit transition of an infinite space with half-spaces 
different in properties substantially and work (lp| imple- 
ments boundary conditions for Wiener processes in path 
integrals. Papers [ill E3 develop a rather sophisticated 
moment technique for tackling the Fokker-Planck equa- 
tion with mixed boundary conditions based on a special 
moment truncation scheme. 

In the present paper we shall extend the method of 
deriving the Fokker-Planck equation from the Chapman- 
Kolmogorov equation in such a way that simultaneously 
consistent boundary conditions can be formulated. Our 
approach is based on introducing physical models for the 
stochastic behavior close to the boundary. We explic- 
itly demonstrate that boundaries break the symmetry of 
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the random forces leading to boundary singularities in 
the Kramers-Moyal expansion. The cancelation of these 
singularities yields the appropriate boundary conditions. 
We explicitly derive the boundary conditions for a re- 
flecting or absorbing barrier as well as boundaries with 
mixed properties, and describe the general procedure for 
the derivation of the boundary conditions for the case of 
the fast diffusion layer. It should be noted that a simi- 
lar anomalous effect of the region boundaries on random 
processes was analyzed in papers [IH 0, [l5[ in numeri- 
cal implementation of Wiener processes in their vicinity. 
Besides, paper [16| applies also to the concept of the sym- 
metry breakdown caused, however, by external fields in 
constructing a generalized master equation for the classic 
and anomalous diffusion processes. 

In principle the present approach can be extended to 
anomalous transport phenomena, e.g., sub- and super- 
diffusion, which are modeled by fractional diffusion oper- 
ators. It is well-known that the formulation of boundary 
conditions for these processes is still a challenging prob- 
lem althou gh several approaches have been developed 
[l7L [H, Il9l l20j] . The procedure outlined in the present 
paper might be helpful in formulating appropriate bound- 
ary conditions for these more involved processes. 

The paper is organized as follows. Section HI1 presents 
the problem under consideration and sketches out deriv- 
ing two types of the Fokker-Planck equations based on 
the general Chapman-Kolmogorov equation for Marko- 
vian processes. Finally it formulates the problem of the 
corresponding boundary conditions and derives the gen- 
eral expressions that should be fulfilled at the boundaries 
of medium. Section IIIII discuses the types of medium 
boundaries and their properties to be taken into account. 
Section IIVI introduces the equivalent lattice description 
of the continuous Markovian process that enables us to 
calculate anomalous kinetic coefficients in the vicinity of 
the boundary. Section [V] is actually the main part of the 
paper, it calculates the boundary singularities. The re- 
sults are used in Sec.|VI]to obtained the desired boundary 
conditions for the forward and backward Fokker-Planck 
equations. 



II. THE CHAPMAN-KOLMOGOROV AND 
FOKKER-PLANCK EQUATIONS 

We consider stochastic dynamics of a Markovian sys- 
tem represented as a point r belonging to a certain do- 
main Q in the Euclidean M-dimensional space R . The 
domain Q is assumed to be bounded by a smooth hyper- 
surface T. When the detailed information about possible 
trajectories {r(i)} of the system motion is of minor im- 
portance the conditional probability called also the Green 
function 

G(r,t\r ,t ) := V{r , t => r, t) 

gives us the complete description of system evolution. By 
definition, the Green function is the probability density 



G(r,t|r ,to) G(r,,t, | r ,t a ) G(r,i|r„,i„) 

> = >b* > 



rOi to 



r, f fa, to r,,t, r.t 

backward FPF ; ■ Forward FPF > 



FIG. 1: Diagram of the Chapman-Kolmogorov equation. The 
symbol l±l denotes summation over the intermediate point r* 
and the arrows illustrate the limit cases t» — > to + and 
t, — > t — matching the backward and forward Fokker-Planck 
equations. 



of finding the system at the point r at time t provided it 
was located at the point ro at the initial time to- 

Since Markovian systems have no memory the Green 
function G(r, £|ro, to) obeys the integral Chapman-Kol- 
mogorov equation that represents transition of the sys- 
tem from the initial point r to the terminal one r within 
the time interval (to, t) as a complex step via an interme- 
diate point r, € Q at a certain fixed moment of time i* 
with succeeding summation over all the possible positions 
of the intermediate point (see, e.g., Ref. [l|) 



G(r,i|r ,to) 



fir* G(r,i|r*,i*) G(r„, i»|r , t ) 



(1) 

The time i» may be chosen arbitrary between the ini- 
tial and terminal time moments, t* € [£()>*]■ Figure [1] 
visualizes this equation. 

Since the domain boundary T is considered to be a 
physical object special properties will be ascribed to it 
and itself can affect the system, for example, trapping it. 
So the symbol of triple integral is used in equation (fTj) 
to underline this feature and where appropriate it should 
be read as 



dr . . . = / dr 



ds . 



ds . 



where the symbol Q + denotes the internal points of the 
domain Q, the boundary T is split from the medium bulk 
because it can differ essentially from the medium bulk in 
properties, and the boundary traps T tr are singled out 
and treated individually by the same reasons. To simplify 
notations a similar rule 



dr. 



dr . 



ds . 



Q+ 



is also adopted. Such a split of integrals is for treating 
motion of the system inside the internal points Q + , its 
possible anomalous transport along the boundary T, and 
the trap effect individually. Besides, according to the 
probability definition, the equality 



drG(r,t|r ,t ) = 1 



(2) 
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holds when the integration runs over all the possible 
states of the system including the boundary traps T tr . 

In what follows a rather general model for the medium 
boundary will be studied. Hear we paid attention only 
to the fact that the boundary traps have to be treated 
individually because the system after been trapped can- 
not leave the boundary reaming in a trap forever. As 
a result if the point i"o belongs to a trap, then for any 
internal point r of the domain Q the Green function is 
equal to zero: 

G(r,*|r 0> to) = for r G T tr , r e Q+. 

Further on the Green function G(r, f|ro, to) for the inter- 
nal initial and terminal points ro, r € Q + will be consid- 
ered. Therefore the general Chapman-Kolmogorov equa- 
tion (TT]) can be reduced by eliminating the integration 
over the traps, so becoming 



A. The backward Fokker-Planck equation 

To implement the limit i» — ► to + let us choose an 
arbitrary small time scale r and consider the Chapman- 
Kolmogorov equation for t* = to+T and an internal point 
ro. Then according to the adopted assumptions the first 
multiplier G(r, t|r», t*) on the right-hand side of (O is a 
smooth function of both the argument r„ and t* whereas 
the second one G(r*, i*|ro, to) exhibits strong variations 
on small spatial scales. So we can expand the function 

G(r,t|r +R, t + r) 

in the Taylor series with respect to the variables r and 
R = r* — ro . The required accuracy is the first order in 
the time step r and the second order in R because the 
characteristic spatial displacement of the system during 
time t is of order r 1 / 2 . Within this accuracy it is 



G(r,t|r ,to)= f[dr.G(T,t\r„,U)G(T.,U\To,ta). (3) G{r, t\r Q + R, t + r) = G(r, t\r , t ) 



+ T- 



In equation ([3]) this elimination is pointed out by the ab- 
sence of one integral matching the traps, cf. the general 
formulation (JTJ) of the Chapman-Kolmogorov equation. 
Within the given integration rule the equality matching 
identity @ is violated, instead, we have 



8G(r f- t ° ) ^ifV?GM|,o,t ) 
dt ^ 



drG(r,t\r ,t ) = 1 



ds tr G(s t r,t\ro,to) < 1 



where the symbol s tr stands for the boundary trap lo- 
cated at the point s 6 T. 

In order to obtain the Fokker-Planck equations two 
additional assumptions must be adopted. The former is 
the short time confinement meaning that on small time 
scales the system cannot jump over long distances or in 
terms of the Green function its first and second moments 
converge and 



1 M 

+ -^Ji'ffV»V»G(r,i|r„,(„), (6) 

where the operator V? = d/dx Q acts only on the argu- 
ment ro of the Green function. The substitution of ex- 
pansion ^ into the Chapman-Kolmogorov equation Q 
reduces it to the following 

9G(r,<|r ,t ) „, , s ni , s 

-t = -9t(r , t , t) G(r, t\r , t ) 



dt 



M 



^il l (ro,io,r)V^G(r,<|r ,to) 



lim 

«->t +o 



drG(r,t\r ,to)\r-r \ p = 0, p=l,2. 



The latter is the medium local homogeneity. In other 
words, the medium where the Markovian process devel- 
ops, i.e. the domain Q should be endowed with char- 
acteristics being actually some smooth fields determined 
inside Q + or at T individually. As a result the Green 
function G(r,t|ro,io) has to be smooth with respect to 
all its arguments for t > to and r, ro G Q + . 

Because the intermediate time t* entering the Chap- 
man-Kolmogorov equation is any fixed value between the 
initial and terminal time moments, to < i* < t, there is 
a freedom to choose it for special reasons. In particular, 
the passage to one of the limits f* — > to + or t* — > 
t — gives rise to either the backward or forward Fokker- 
Planck equation, respectively (Fig. [1} . 



M 

+ £ y '(ro,to,T)V°V J G(r,t|r ,io), 

(7) 

where the quantities 

9t(ro > t ,r) = l- // rfRG(r + R, t + r|r , t ) , (8) 



il l (r ,t ,r) = JJ dRR l G(r + R,t + T\r ,to), (9) 

Q 

£ y '(r ,t ,r) = 1 J J dRfl i JJ*G(r + R J *o+T|ro,to) 

Q 

(10) 

have been introduced. Besides, the first term on the 
right-hand side of {7} has been assumed to be small and 
tend to zero as r — > which is justified based on the 
results to be obtained. 
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For an internal point ro and, thus, separated from the 
boundary T by finite distance the time step r can be cho- 
sen so small that it is possible to construct a neighbor- 
hood of the point ro with the following properties. First, 
deviation of the Green function G(ro + R, to + r|ro,to) 
from zero outside this neighborhood is ignorable due the 
first assumption about the short time confinement. Sec- 
ond, inside it the medium can be regarded as the homo- 
geneous space K M by virtue of the second assumption on 
the local homogeneity. In this case actually replicating 
the proof of the Law of Large Numbers using the gen- 
eration function notion (see, e.g., Ref. it is possible 
to demonstrate that quantities (HJ) and (fT0|) scale linearly 
with t. The difference of quantity ([8]) from zero is ignor- 
able. Therefore for internal points we can introduce the 
drift velocity z/(r, t) and the diffusion tensor D %3 ; (r, t) by 
the expressions 



symmetry breakdown of 
random walks near 
the medium boundary 
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FIG. 2: The effect of the boundary impermeability on the 
Markovian system motion. Schematic illustration. 



v*(r,t)= lim - / dRB l G(r + R,t + T\r,t), (11) 

T J 

Q+ 

D ij (r,t) = lim — / dRR l R j G{r + R, t + rlr, t). 
t^+o 2r J 

Q+ 

(12) 

Then for the internal points the division of equation (JjJ 
by t and the succeeding passage to the limit r — > +0 
yield the backward Fokker-Planck equation 

dG(r, t|r ,t ) - , , 
^7 = £fp b {&(r,t|r ,to)j , (13) 

where the backward Fokker-Planck operator is 

M M 

C FPb := ]T ^'(ro,to,r)V°V J + ^i; i (r ,to,T)V°. 

(14) 

We note that the backward Fokker-Planck equation acts 
on the second spatial argument of the Green function 
G(r,t\r ,t ). 

This Fokker-Planck equation should be supplemented 
with the initial condition and the boundary condition. 
By construction, at the initial time to the system was 
located at the internal point ro, so the initial condition 
just writes the Green function in the form of the Dirac 
(^-function 

G(r,t\r ,t )\ t=to =6(r-r ), (15) 

The boundary condition interrelates the values of the 
Green function and its derivatives at the internal points 
adjacent to the domain boundary T, i.e. values ob- 
tained by continuation ro — > s from some internal point 
ro G Q M+ to a boundary point s e T. 



B. The boundary condition problem for the 
backward Fokker-Planck equation and the vector of 
boundary singularity 

The direct implementation of the passage to the 
boundary points, however, raises a certain problem. Ex- 
pansion ([7]) exhibits irregular behavior within the joint 
passage to limits r — > +0 and r — > s. When the for- 
mer r — > +0 precedes the latter ro — > s no boundary 
conditions are got at all. 

In the opposite order, i.e. when the passage ro — > s is 
performed first, the kinetic coefficients (|5j)- (fTU|) change 
the scaling type; now they vary with time r as yfr at the 
leading order. The matter is that a path of Markovian 
system is not smooth at every point and its characteristic 
variations on small time scales about t are proportional 
to s/t. For the internal points of the domain Q the path 
deviations in opposite directions are equiprobable within 
accuracy y/r. As a result the coefficient il*(r, t, r) be- 
comes a linear function of the argument r. In some sense 
the given anomaly in the Markovian dynamics is hidden 
at the internal points and reflected only in the linear r- 
dependence of the second order moments £ y (ro, to 7 T ) of 
the Green function G(r, t|ro, to). The medium bound- 
ary T breaks down this symmetry because, in particu- 
lar, it prevents the system from getting the points on 
the opposite side. Since the system displacement re- 
mains the same magnitude the terms il z (r, t,r) acquire 
the root square dependence on the argument r. In a cer- 
tain seance the medium boundary reveals this anomaly 
(Fig. The succeeding division of expansion ([7]) by r 
gives rise to singularities of the type r -1 / 2 which will be 
referred to as boundary singularities. 

The medium boundary can affect the system dynam- 
ics in a more complex way, here, however, we currently 
confine our speculations only to the effect of its imperme- 
ability Since the boundary T confines the system motion 
only in the normal direction it is quite natural to expect 
that the boundary singularities quantified in terms of di- 
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verging components II 1 (s, t, t)/t will form a vector object 
b that is determined by mutual effect of two factors. The 
former is the spatial orientation of the medium bound- 
ary T described by its unit normal n. The latter is the 
spatial arrangement and intensity of random Langevin 
forces governing stochastic motion of the given system. 
They are characterized by the diffusion tensor D % i (r, t) . 
Within a scalar cofactor we have only one possibility to 
construct the vector b = {b 1 } using the two objects, 

M 

3=1 

or, in the vector form 

b = D n. (17) 

The validity of this construction will be justified in the 
present paper and b will be referred to as the vector of 
boundary singularities. To be rigorous it should be noted 
that in the general case the correct expression for the vec- 
tor of boundary singularities should use the operator 
obtained from the diffusion tensor D 1 ^ by lowering one of 
its indices, namely, b l = J2j D]ni (such details are dis- 
cussed in Sec. IIV A[) . However dealing with orthonormal 
bases as it is the case at the initial stage of the current 
consideration the tensors D 1 - 7 and D* coincide with other 
in the component magnitudes. So in order not to over- 
load the reader perception and the mathematical con- 
structions expressions similar to (| 16[) will be used where 
appropriate. 

The notion of the boundary singularity vector enables 
us to write immediately the desired boundary condition 
when the medium boundary just confines the system 
motion. In this case the first and third terms on the 
right-hand side of expansion ([7]) are absent and the cor- 
responding singularity caused by the sequence of transi- 
tions ro — > s and then r — ► takes the form 

1 M 

-=5>*(s)V° G(r,*|ri )| ro _ 

V i — 1 

i M 

= ~r E ^'(MoK(s)V? G(r,i|r ,t )| ro ^ s • 

Naturally, for the internal points r and ro the Green func- 
tion G(r, i|ro,io) cannot exhibit any singularity whence 
it follows that the cofactor of the singularity t -1 / 2 must 
be equal to zero, i.e. 

M 

J2 ^'MoK(s)V° G(r,t|r o ,i )| ro ^ s = 0. 

i,3*=l 

It is the very know expression for the boundary condi- 
tion of the backward Fokker-Planck equation which is 
typically obtained in another way applying to physical 
meaning of the Green function (see, e.g., Ref. [l[). 



The present paper is devoted to deriving the boundary 
conditions for the Fokker-Planck equation applying to 
the notion of the boundary singularities. A more general 
situation will be studied justifying also these qualitative 
speculations. Currently we can state that the bound- 
ary condition for the backward Fokker-Planck equation 
should stem from the requirement for the boundary sin- 
gularity terms to vanish in expansion ([7]), i.e. when 
r <£ T 

-*m(r ,t Q ,T)G(r,t\r ,t ) 

M 

+ ^*if(ro,<o,T)V°G(r,i|r ,i ) 

M 

+ J2 *£ ij (ro,t ,T)V° i V° j G(r,t\r ,to)=0, (18) 

where the symbol * labels the components of the corre- 
sponding kinetic coefficients scaling as r 1 / 2 . It should 
be pointed out that in expression (I18|) the argument ro 
is an arbitrary point of a thing layer T T adjacent to the 
boundary T, which is designated by the symbol <s. When 
t — * its thickness also tends to zero (as t 1 / 2 ), however, 
before passing to the limit r — ► the layer T T remains 
volumetric. 

Now let us discuss similar problems with respect to 
the forward Fokker-Planck equation matching the other 
possibility of passage to the limit case in the Chapman- 
Kolmogorov equation ([3]). 



C. The forward Fokker-Planck equation 

The Chapman-Kolmogorov equation © also allows for 
the limit where the intermediate point tends to the termi- 
nal one, i.e. i* = t — r with r — * +0. In this case the for- 
mer cofactor G(r,t|r*,t — r) on the right-hand side of Q 
exhibits strong variations on small spatial scales whereas 
the latter one G(r*,t — r|ro, to) becomes a smooth func- 
tion of the argument r*. Now, however, applying directly 
to an expansion similar to that have been used in deriving 
the backward Fokker-Planck equation is not appropriate. 
The matter is that in this way the integration runs over 
the initial point r* of the Green function G(r, t|r*, t — r) 
and appearing coefficients similar to quantities (15j)- (fTU|) 
have another meaning. In particular, an integral similar 
to Q can deviate from unity essentially. 

To overcome this problem the Pontryagin technique is 
applied [2lJ. It is rather similar to the Kramers-Moyal 
approach (see, e.g., @) but is more suitable for tackling 
the boundary singularity. Let us consider at the first 
step some arbitrary smooth function </>(r) determined in 
the domain Q and integrate with it both the sides of the 



6 



Chapman-Kolmogorov equation ([3]). In this way we get 



J j dr^(r)G(r,4 +r|r ,t ) 



drdr* 4>{y) G(r,t* + r|r*,i*) G(r*,t*|ro,*o) ■ 

(19) 

For a rather small time scale r the Green function 

G(r,t* +r\r*,t*) 

is practically located within some small neighborhood of 
the point r* . Thereby the function </>(r) can be expanded 
in the Taylor series near the point r* with respect to the 
variable R = r 

M M 

0(r)-0(r,) + ^^V^(r,) + - ]T i^V*V*0(r*) . 

Beside, since the Green function G(r,t >t + r|ro,io) de- 
pends smoothly on r the expansion 



G(r,i* +r|r ,i ) = G(r, i*|r , i ) + t 



is also justified for a small value of r. 

Then the substitution of the last two expressions into 
equation (|19p with succeeding integration over R and the 
replacement of the dummy variable r* by r as well as t» 
by t yields 



drc[>(r) 



_dG{r,t\r ,t ) 



Of 



dr{ <j>(r) 



M 



■ J2 v * v ^( r ) 



-«R(r,t,T)G(r,t|r ,to) 
il J (r,i,r)G(r,i|r ,t ) 
£«(r,t ) r)G(r,t|r ,to) 



(20) 



Here the coefficients 9i(r,i, r), iT(r, £,t) and £ y (r, t,r) 
again exhibit anomalous behavior within a narrow layer 
T r adjacent to the medium boundary T (Fig. [3J. As 
should be expected and in accordance with results to be 
obtained the thickness of this layer scales with time r 
as t 1 / 2 . These coefficients themselves also scale as r 1 / 2 . 
As a result the corresponding part of integral (|20|) scales 
as r. Thereby after dividing both of the sides of (|2"01) 
by r with the following passage to the limit r — > the 
contribution to (|20|) caused by integration over this layer 
remains finite. Therefore to analyze the properties of the 
integral relation (12"0)) the domain Q is split into this layer 



layer of 
boundary \Tf r 
singularitics 




FIG. 3: Structure of integral (|20|l and division of the region 
Q into the layer of boundary singularities and the internal 
points with regular behavior of the kinetic coefficients. 



of boundary singularities and the internal part. After the 
passage to the limit r — ► this division matches treating 
individually the boundary T and the internal points Q + . 

Keeping the aforementioned in mind the integral ex- 
pression (|20p is represented as a sum of two terms, the 
integral over the layer T r denoted with the formal sym- 
bol of surface integral and the integral over the internal 
part Q + of the domain Q 



// 



dr . . . = d> dr 



dr. 



(21) 



Q+ 



Let us consider the second term first. Inside the region 
Q+ the kinetic coefficients iP(r, t,r) and £ y (r, t, t) be- 
have in regular way, i.e. they scale as r according for- 
mulae (fTTj) and (fT2")) . whereas the term 9t(r, t, r) vanishes 
at all. So dividing the corresponding part of the integral 
relation ([2H)l by r and passing to the limit r^Owe have 



drcf>(r) 



dG{r,t\v ,t ) 



Of 



J drj^V^(r) 



A I 



v t (r,t,T)G(r,t\r ,t ) 



D ij (r,t, T )G(r,t\r ,t ) 



(22) 



Using the Gauss divergence theorem this integral in turn 
is split into two parts, surface and volume ones: 



dr . . . = d> ds . 



dr. 



(23) 
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The volume integral has the form 
"5G(r,*|r ,to)" 



integral over the layer T T of boundary singularities 



dr<fi(r) 



dt 



J dr0(r)j -^V 



M 



w l (r,t,r)G(r,i|r ,t ) 
D«(r s t ) r)G(r,t|r 0) to) 



(24) 



The latter equality immediately gives rise to the forward 
Fokker-Planck equation. 

Indeed, currently <f>(r) is an arbitrary smooth function 
and no addition constrains will be imposed further on it 
for the internal points of the domain Q. So applying to 
local variations of (f>(r) at an arbitrary internal point r 
(Fig. [3J we see that the left and right sides of flM]) should 
be equal to each other for the points r <= Q + individually, 
getting the forward Fokker-Planck equation 



dG{r,t\r ,t ) 
dt 



= £ FPf {G{r,t\r ,t )} 



(25) 



with the forward Fokker-Planck operator 

M 



A I 



i=i 



(26) 



Here the symbol stands for a function acted by this 
operator. It should be also pointed out that the Fokker- 
Planck operator acts on the first spatial argument of the 
Green function. 

The forward Fokker-Planck equation can be also writ- 
ten in the conservation form 



dG(T,t\r ,t ) 
dt 



M 



^V i J i {G(r,t|r 0) to)} = 0, (27) 



A I 



j> dr G(s, t\r ,t ) J jr Vi^(s) *if(r, t, r) 

Y I. i— 1 

M \ 

- cf>(s) *<K(r, t,r)+J2 ViV^(s) *£ y (r, i, r) . (29) 
i,j=i J 

Here the symbol dr as well as presence of the argument 
r in the singular components of the kinetic coefficients 
takes into account the fact that before the passage to 
the limit r — ► the layer T r is volumetric. The Green 
function G(r, t|ro,to) as well as the test function cf>(r) 
and its derivatives exhibits minor variations across the 
layer T T so their argument r have been replaced by the 
corresponding nearest point s laying on the boundary T. 

The latter term is due to the part of expression (|2"2")l re- 
maining after integration using the convergence theorem 
and can be written in the form 



/M 
ds V J 0(s)n i (s) 



£>y(s > f > 7-)G(s J i|ro,to) 



. M 

= -& ds0(s) 2 n\s) r {G(s, t|r , to)} , (30) 

j. 2—1 

where n(s) = {n*(s)} is the unit normal to the boundary 
T at point s directed inwards the domain Q. 

Leaping ahead we note that the appropriate choice of 
the boundary values of the test function </>(s) and its 
derivatives fulfils equality (|29p and, at the next step, gives 
rise to the required boundary condition for the forward 
Fokker-Planck equation. Let us demonstrate this for the 
impermeable boundary using the notion of the boundary 
singularity vector b. Namely, we again assume that for 
an internal point r located in the vicinity of a boundary 
point s, i.e. r <s s 

lT(r,i,T) ocb l (s) =J2D ij (s,t)u j {s). 



with the probability flux operator J = {J 1 }*^ 



J l {0} 



7 = 1 ^ ' 



+ </(r,;,r)0. (28) 



The forward Fokker-Planck is naturally supplemented 
with the same initial condition (fl~5l). 



D. Boundary relations for the forward 
Fokker-Planck equation 

Splits (|2"Tj) and (|2"5)) give rise to two additional terms. 
The former one is related to the first split and is the 



In this case only the first term in equality (|29[) remains 
and it is fulfilled when 



^D«(s,ty(sMs)=0. 

ij=l 



(31) 



Equality (|3Tj) just relates the boundary values of the test 
function cf>(&) with its derivative along the boundary nor- 
mal n(s). So for an arbitrary smooth function c/f>T(s) 
determined at the boundary T it is possible to construct 
the appropriate function </>(r) determined in the domain 
Q and meeting equality (pH]) (see Fig. [3]) . So in the given 
case the left-hand side and, thus, the right-hand side of 
expression (1301) becomes zero. Since the integral on the 
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impermeable 
boundary 



absorbing 
boundary 



boundary with 
fast diffusion layer 




The third type boundaries are widely met, for example, 
in polycrystals or nanoparticle agglomerates. The grain 
boundaries contain a huge amount of defects and as a re- 
sult the diffusion coefficient inside the grain boundaries 
can exceed its value in the crystal bulk by many orders. 
Therefore impurity propagation in polycrystals is gov- 
erned mainly by g rain boundary diffusion (for a review 
see, e.g., Ref. 22[ and references therein). In terms of 
random walks the effect of the fast diffusion layer is re- 
duced to extremely long spatial jumps made by an walker 
inside it. It is natural to characterize such a boundary 
layers by its thickness A about the atomic spacing and 
the ratio of the diffusion coefficients inside the boundary 
layer and in the regular crystal lattice q ^> 1 . 



FIG. 4: Three types of boundaries under consideration. 

right-hand side of (|30[) contains an arbitrary function 
4>(s) determined at the boundary T the equality 

M 

5> l (s)T {G(s,i|r o ,M} = (32) 

i=l 

holds for every point of the boundary T individually. 
This expression meaning the zero value of the probability 
flux in the direction normal to the boundary T matches 
well the physical seance of its impermeability. 

However, to derive the boundary conditions for the 
Fokker-Planck equations more sophisticated construc- 
tions are necessary. Besides, in order to take into ac- 
count other possible properties of the medium boundary 
its model should be specified. 

III. BOUNDARY TYPES 

In the present paper, to be specific, we consider three 
typical examples of medium boundaries. They are (?) the 
impermeable boundary, [n) the boundary absorbing par- 
ticles, and (m) the boundary with a thin adjacent layer 
characterized by extremely high values of the kinetic co- 
efficients, the fast diffusion boundary (Fig. |4j). 

The first type matches a medium whose boundary is 
similar to its bulk in properties, the boundary points dif- 
fer from internal ones only by the absence of medium 
points on one side. As a result a random walker hopping 
over the medium points just cannot pass through the 
boundary returning to the medium bulk after getting it. 

The second type is similar to the first one except for 
the fact that the walker can be trapped at the boundary 
and will not return to the medium anymore. In this case 
the corresponding boundary conditions are typically used 
in describing the first passage time problem or diffusion 
in solids with fixed boundary values of impurity concen- 
tration C s (see, e.g., Ref. [l|). Generally the boundary 
absorption is described by the rate aC s , where a is a 
certain kinetic coefficient. 



IV. EQUIVALENT LATTICE 
REPRESENTATION OF RANDOM WALKS 
NEAR THE MEDIUM BOUNDARY 

The derivation of the Fokker-Planck equations, the for- 
ward and backward ones, requires calculation of three 
quantities 5H(r, t,r), iP(r, t, r), and £ w (r, t,r) specified 
by expressions ©-(HUJ). They are the moments of the 
system displacement R during the time r treated as an 
arbitrary small value. In order to obtain the desired 
boundary conditions these quantities should be found in 
the vicinity of the medium boundary T or, more pre- 
cisely, in its neighborhood T T of thickness about {Dt) 1 / 2 , 
where D is the characteristic value of the diffusion tensor 
components. To study the boundary effects it suffices to 
consider a rather small region wherein the medium and 
its boundary are practically homogeneous in properties 
and, in addition, the boundary geometry is approximated 
well by some hyperplane. In this region the system mo- 
tion will be imitated by random walks on a lattice con- 
structed as follows. 

First, the elementary steps of the random walks on it 
are characterized by a time r such that 

T a « T (33) 

and the arrangement of the lattice nodes, i.e., their spac- 
ings {a-i} and the spatial orientation should give us again 
the same diffusion tensor D as well as the drift field v for 
the internal points on time scales r a « t C t. The indi- 
vidual hops of a random walker between the neighboring 
nodes actually represent a collection of mutually inde- 
pendent Langevin forces governing the random system 
motion in the given continuum. Second, the boundary T 
is represented as a layer of nodes To between which the 
walker can migrate via elementary hops. In other words, 
the aforementioned collection of mutually independent 
Langevin forces has to contain components acting along 
the boundary T and one component moving the walker 
towards or from T. Other characteristics of this effective 
lattice may be chosen for the sake of convenience. At 
final stage we should pass to the limit r a — > returning 
to the continuous description. 
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A. Diffusion tensor representations 

In order to construct the required lattice let us con- 
sider Markovian random walks {r(i)} in M-dimensional 
Euclidean half-space R M+ made of vectors 

r = {x x ,x 2 ,...,x M } 

such that 

M 

r • n :— x z n l > , 

i=l 

where n — {n 1 , n 2 , . . . , n M } is a certain unit vector. The 
boundary of R M+ , i.e. the hyperplane T = {r • n = 0} 
perpendicular to the vector n is, in its turn, the Eu- 
clidean space R M_1 of dimension M — 1. The half-space 
R M+ and, correspondingly, the hyperplane T are as- 
sumed to be homogeneous. The latter means the local 
properties of the random walks under consideration to be 
independent of position in space; naturally the boundary 
and internal points are not equivalent. In particular, the 
diffusion tensor D and drift vector v are the same at all 
the internal points of the half-space R M+ . 

In this case the components of the drift vector and 
diffusion tensor are determined by the expressions (cf. 
formulae (j8l)- (fT0| ) 

v* = ±(5X*(t,T)), (34) 

D ij [SX% t) - v*t] [8X 3 {t, t) - v j r] J> . (35) 

Here the random variable SX l (t,T) := x' l (t + r) — x^if) 
and r = {x 1 } is an arbitrary internal point, the ob- 
servation time interval r should be chosen to be small 
enough that the length scale (Dt) 1 ^ 2 be much less then 
the distance between the point r and the boundary T, i.e. 
Dt <C (r • n) 2 , and the triangular brackets (. . .) stands 
for averaging over all the random trajectories passing 
through the point r at time t. It should be noted that 
due to the space homogeneity the passage to the limit 
t — ► can be omitted which is necessary in the general 
case. 

In what follows nonorthogonal bases will be used. So, 
keeping in mind the tensor notation (see, e.g. Ref. [HI), 
the upper and lower indices will be distinguished. In 
these terms {x 1 } or just x % is a vector, whereas, the col- 
lection of the basis vectors e$ is a covector. According 
to definitions ([33)) and l|34p the objects D 1 ^ and v l are 
contravariant tensors. In addition, if the basis e has the 
form e = ex © e, where ex is the basis of the hyperplane 
T and the vector e does not lay in it, then the Greek 
letters will label the tensor indices corresponding to the 
hyperplane T to simplify perceiving this fact. 

In order to deal with the diffusion tensor in a 
nonorthogonal basis e = {e^} the metric tensor is also 
necessary. It is defined as 

9ij •■= (?i • *j) (36) 



and is the kernel of the scalar product of two vectors r 
and r, namely, 

M 

(r • r) := ^ g ijX l x ] . 

For an orthonormal basis the metric tensor g^j = 5ij, 
where 5y is the Kronecker delta. The metric tensor g^j 
defines the conversion of contravariant tensors into co- 
variant ones, in particular, 

M M 

D)=J2 Dik 9kj, Di = Y,9ikD k \ (37) 

k=l k=l 

as well as 

M 

= J2 9ik 9j P D kp . (38) 

k,p— 1 

Due to the diffusion tensor D IJ as well as the metric ten- 
sor gij being symmetric the tensor Dij is also symmetric, 
whereas the tensors £J l ■ and are identical and so de- 
noted further as £)*• . The tensor can be regarded as a 

certain operator T> acting in the space R and the tensor 
Dij specifies a quadratic form 

M M 

r-Vr= I H-' J = D l3 x l x i > . (39) 

ij,k—l ij — 1 

The quadratic form (|39[) is positive definite. To demon- 
strate this a random variable 

M 

SL = Y^ [SX p {t, t) - v p t] (e p -£) 
P =i 

M 

= J2 [SX p (t,r)-v p r]g pi e 

p,i—l 

is considered, where £ — YlfLi e i^ 1 1S an arbitrary vector 
in the space R M and the metric tensor definition |33)) 
has been taken into account. Whence we have a chain of 
equalities 

M 

0<([<5L] 2 ) = 9 Pt g P H>i l i 1 ' 
x ([5X p (t,T) -v p t] [SX p '(t,T) -v p 't]} 

M 

= Y, 2tL>pp ' dpigp'fir' 

p^p' — 1 
M M 

= Y IrDadH 1 ' = Y ItD^'UUi. 

i.i' — l p.p' — 1 
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So for any arbitrary vector l l and covector Zj the inequal- 
ities 



M 



M 



E DijPp > o, 



E l) ' i, -'> 



> 



(40) 



hold. The covector and vector representations of the 
same object are related as Z^ = y^-_i ggP; within or- 
thonormal bases they are identical. 

Due to the symmetry of the tensor Dij and the quadra- 
tic form (|39|) being positive definite all the eigenvalues of 
the operator V are real positive quantities and its eigen- 
vectors form a basis in the space R M which can be chosen 
to be orthonormal one, see, e.g., Ref. [24[. In this basis 
the diffusion tensor takes the diagonal form. Thereby 
the corresponding eigenvectors and eigenvalues specify 
the directions and intensity of the mutually independent 
Langevin forces governing random walks in the medium 
under consideration. Unfortunately, in the general case 
where all the eigenvalues are nondegenerate this basis is 
unique. So it cannot be used in constructing the desired 
lattice in the vicinity of the medium boundary T because 
one can meet a situation when none of the basis vectors 
is parallel to the hyperplane T. In order to overcome this 
problem we will construct a special nonorthogonal basis 
applying to the following statement. 

Proposition 1 Let M. M+ = {r-n > 0} be a homogeneous 
half-space bounded by the hyperplane T = {r-n = 0} and 
e = {ei,e2, . . .enj} be a fixed arbitrary basis o/R M . In 
this basis the components of the diffusion tensor 
as well as the matric tensor {gij} are given. Then there 
is a basis b — bx © bu with the following properties. 

First, it is composed of a certain orthonormal basis bx 
of the hyperplane T and a unit vector b&i not belonging 
to T that is determined by the expression 



1 



M 

E ^ D > 3 



(41) 



Here according to the construction of the half-space 
n = {n 1 , n 2 , . . . n M } is the unit vector normal to the hy- 
perplane T and the normalization factor 



to 



M 



i,j,p,k— 1 



1/2 



(42) 



Second, in the basis b the diffusion tensor takes the diag- 
onal form 



\D l 



V 1 
V 2 





V 



M 



(43) 



where all its diagonal components are positive quantities, 
{T>i > 0}, with the value T>m being given by the expres- 
sion 

-l 



M 



E D ^ 



(44) 



Third, let, in addition, the initial basis be of the form 
e = ex © n, where ex = { e a}i /_1 is a certain basis of 
the hyperplane T, and Ut = \\u a a\\ be the transformation 
of the hyperplane T mapping the basis bx onto the basis 

ex, i.e., bx ► ex- By mapping | — > n the transfor- 
mation Ut is complemented to a certain transformation 
U of the compleat space R M , namely, if r is an arbitrary 
vector of the space K M with the coordinates specified by 
its expansion over the bases e and b: 

Af-l ftf-1 

r = e 7^ 7 + nxM = E h ^ + bM< ^ M > 

7—1 7=1 

then its coordinates are related by the expressions 



M-l 

E' 

7=1 

c M = w 



D 



MM 



D1 M, 



.At 



JjMM ' 

and for the inverse transformation 

1 



M-l 

E-% 

7=1 
j~)MM 



-D 



aM i-M 



■ M 



(45) 

(46a) 
(46b) 

(46c) 
(46d) 



Here U„ 



\u a a\\ is the operator inverse to the oper- 



ator Z7 T , i.e. meeting the identity 
Besides, the equality 



M-l 

l 



U 7 U P ~ °P- 



M-l 

E U^U^ 
7=1 



= D a0 



1 



D 



MM 



jjaM jjpM ^ 



holds. 



Since the proof of this proposition requires just formal 
mathematical manipulations related weakly to the sub- 
ject matter of the paper it is presented in the individual 
Appendix [5] 

Comments on Proposition [I] First, it is worthwhile to 
note that the basis vector bM constructed by expres- 
sion (|4"Tj) is actually the vector b of boundary singularities 
(expression (|16p ) normalized to unity. 

Second, for the initial basis e of the general form actu- 
ally expression (147[) persuades us to introduce the surface 
diffusion tensor 



(48) 
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that describes the system random motion along the hy- 
perplane T. Indeed, in a basis bx ® n the components of 
this tensor belonging to the hyperplane T coincide with 
ones given by expression (|47[) and are equal to zero when 
one of its indices matches the vector n. 

Third, when the initial basis e is orthonormal the ex- 
pressions of Proposition [T] can be simplified. Indeed, in 
this case the matric tensor := (e$ -e,-) = 5jj is the unit 
matrix and it is possible not to distinguish between the 
upper and lower tensor indices, in particular, all the com- 



ponents Di 



n 7 ~J 
the initial basis has the form e 



D\ = D lJ are identical. If, in addition, 
ex© n expressions (j4l~j) - 



44l) become 



>M 



M-l 

^ e 7 D 7 A/ + tiDmm 

.7=1 



(49) 




where the coefficient 



M-l 
L7=l 



-,1/2 



lM + D MM 



(50) 



and the value 



M-l 



v M = 



D 



— E d *m 

MM 



D 



MM ■ 



(51) 



7=1 



Besides, the inverse transformation matrix ||u a/ g|| coin- 
cides with the direct transformation matrix transposed, 
i.e. u al3 = up a . □ 

Proposition [1] prompts us to use the basis b = {hi} 
in describing random walks in the half-space R M+ . For 
its internal points the continuous random walks are rep- 
resented as a collection of mutually independent one- 
dimensional Markovian processes {C(t)} 



r(t) = biCit) 



b, / dt' £*(*'), 
'o 



(52) 



where the Langevin random forces {C{t)} meet the cor- 
relations 



(e(t))=v\ (53) 
?(t)e'(t'))=2V i 6u'S(t-t'), (54) 



and {v 1 } are the components of the drift velocity v = 
hiV 1 in the basis b. As could be shown directly these 
random forces lead to expressions (fM)) and (|3"S"|) . 



B. Equivalent lattice random walks 

The desired lattice is constructed as follows (see also 
Fig. [5] for illustration) . At the first step a set of nodes 



FIG. 5: The lattice random walks imitating continuous 
Markovian process in the half-space R 3+ . Here T is the 
boundary of R 3+ , the axes x 1 , x 2 are chosen to be directed 
along the vectors bi, D2 of the basis bx, the axis x 3 is normal 
to the plane T, whereas the basic vector 63 is not normal to 
it in the general case. The values 01,02,03 are the lattice 
spacings and grey arrows show possible hopes to the nearest 
neighbors. 



{ax} is fixed on the boundary T such that 



ar(nr) = h a a, 



M-l 

E 

a=l 



(55) 



where nx = {ni, tt-2, ■ ■ ■ , tim— 1} is a collection of integers 
taking values in Z and the lattice spacings a a are chosen 
to be equal to 



^2T a MV a 



(56a) 



Here r a is any small time scale meeting inequality (|33[) 
and being the time step of lattice random works; an 
walker hops to one of the nearest neighbors in time r a . 
Such jumps are illustrated by grey arrows in Fig. [5] These 
nodes are regarded as the boundary layer To of the lat- 
tice to be constructed. Then the layer To as a whole 
is shifted inwards the region R M+ by the vector Hm^m, 
where 



y/2r a MV. 



M 



(56b) 



Then this new layer Ti in turn is shifted by the same 
vector M M+ , giving rise to the next layer T2 of nodes 
and so on. In this way we construct the system of layers 
{Tfc} making up the desired lattice and exactly random 
walks on this lattice will imitate the continuous process 
in the half-space R M+ . 

Let us now specify the probability of hops from an 
internal node n to one of its nearest neighbors n' along 
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a basis vector by the expression 



1 

2M 



2a 



■v'Xi ■ 



(57) 



Here the random value Xi = il takes into accounts the 
possibility of jumps along the vector b; or in the opposite 
direction. The sequence of such hops with time step r a 
represents equivalently the continuous process rather far 
from the boundary T. Indeed, due to the law of large 
numbers (see, e.g. (25[) two Markovian processes are 
identical if on a rather small time scale both of them 
lead to the same mean and mean-square values of the 
system displacement. By virtue of (157)1 one hop of the 
walker is characterized by the following mean values of 
its displacement Sr = hi5( l 



n' 



(58) 
(59) 



where the sums run over all the nearest neighbors n' of 
the node n. According to and (|35p actually the 
same mean values of the system displacement during the 
time interval T a are given by the continuous random pro- 
cess. Rigorously speaking, the latter mean value and one 
corresponding to the continuous random process are not 
identical, but their difference 

(bi.b^VT^ 

is of the second order in the time scale r a whereas the 
leading terms are of the first order. Thereby choosing 
the time scale r a to be sufficiently small we can make 
this difference ignorable. 



C. Properties of the boundary layer To 

In order to describe the boundary effects on random 
walks special properties should be ascribed to the nodes 
of the boundary layer T . It is worthwhile to noted that 
it is the place where the model of the medium boundary 
does appear for the first time. 

Keeping in mind the boundary types discussed in 
Sec. IIII1 first, each boundary node is regarded as a unit 
of two elements, the lattice node itself and a trap. If a 
walker jumps to a trap it will never return to the lattice 
nodes. The introduction of traps mimics the absorption 
effect of medium boundaries. Second, possible fast diffu- 
sion inside a thin layer adjacent to medium boundaries 
is imitated in terms of multiple steps over the boundary 
nodes during the time interval r . These constructions 
are illustrated in Fig. [5J 

For the walker located at a certain boundary node the 
probabilities of hopping to the internal neighboring node, 
Pi, or being trapped, P tr are specified as 



Pi = 



1 



M 



M ' 



(60) 



T 

/! 

node trap 



T, 



.0 




FIG. 6: Characteristic properties of random walks in the 
boundary layer To. The left inset visualizes possible hops 
from the boundary layer. The main fragment illustrates the 
walker jumps inside the boundary layer To which can be com- 
plex and comprise many elementary hops. The latter feature 
imitates possible fast diffusion inside a certain thin layer ad- 
jacent to crystal boundaries 



where the coefficient a a quantifies the trapping (absorp- 
tion) effect. Leaping ahead, we note that the coefficient 
<7 Q can be assumed to be a small value because its mag- 
nitude a a — > as r a — > within the collection of lat- 
tices leading to the equivalent description of the random 
walks on time scales r a -C t <C r. These probabilities 
have been chosen to constitute the probability of walker 
motion along the direction of the basis vector bja equal 
to the same value for the internal points, 



Pi + Pu 



1 

M 



Therefore the probability for the walker being initially 
at a boundary node nr to make a jump within the bound- 
ary layer T is 



Ft 



M -1 

M 



(61) 



At first, let us consider the case where such jumps are 
the elementary hops to one n T of the nearest neighboring 
nodes in To- Then, following actually construction (15T 
its conditional probability is written as 



P 



(i) 



1 



M 



2(M - 1) (M - 1) 2a a 



VrXa ■ 



(62) 



Here, as before, the value \a = ±1 is ascribed to the 
walker hop along the basis vector h a or in the opposite 
direction, v^- are the components of the drift velocity 
inside the boundary layer in the basis bx- It should be 
noted that regular drift inside the boundary layer and 
the medium can be different in nature, which is allowed 
for by the index T at the boundary components of the 
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drift velocity. The adopted expression ([62]) , as it must, 
obeys equalities similar to expressions (1551) . (f59")l . namely, 
for the displacement <5r T = h a SC a along the boundary T 



V. BOUNDARY SINGULARITIES 



mET 



PT Y P (1) SC rr, SC = 2r a V, 
1 / j ii y in ^n T m ^n T in a 



aO a (3 • 



(63) 



(64) 



The fast diffusion inside the boundary layer T is imitated 
by complex jumps made up of g successive elementary 
hops within the time r a . In this case the walker can get 
not only the nearest neighboring nodes but also relatively 
distant ones. The conditional probability of such a 5-fold 
jump from node n T to node n' T is given by the expression 



P 



(9) _ 



E 



p 



(i) 



X P 



(1) 



mi ,iTi2 ,...,m g _i GT 



•••xP (1) x P y 

in , 111 , 111 , ri x 

g — 2 g — 1 g — 1 T 



(65) 



By virtue of (f6"3"| and (j6"4"|) the probability function P^f ' n , 
of g-fold jumps gives the following values for the first 
and second moments of the walker displacement <5r T = 

J2a=i h a$( a in the la y er T o 



Pr Y P {9) SC r 

m£T 

Pr Y P {9) SC m SC 

1 / < n T m ^n T m ^n T r 



gT a vr , (66) 

2{gT a )V a 5 a i3 
M 



(M - 1) 



a 



(67) 



In expression (|67[) we again have ignored terms of or- 
der gr^ because the displacement of a walker along the 
boundary T caused by its migration inside the layer To 
is considerable only for g 3> 1 as will be seen further. In 
latter case the conditional probability (163)) of transition 



from the node n T to the node 



rix 



({m a } are integers) 



can be approximated by the Gaussian distribution 
M-l 



P 



>(<?) ; = / 



x exp 



M-l 
2ng 



(M - 1) 
25 



M-l 
a=l 



gM 



(M-l) a a 



(68) 



by virtue of the law of large numbers and expres- 
sions (pijl . (JMI), 

The desired lattice random walks imitating the con- 
tinuous Markovian process in the vicinity of the medium 
boundary T is constructed. 



As discussed in Sec. Ill Bl the medium boundary T 
breaks down the symmetry of random walks in its vicin- 
ity, which is reflected in the anomalous behavior of the 
means quantities (|8|)- lfT0|) near the boundary T. To 
quantify this effect it is necessary to calculate the given 
integrals near the boundary T for any small time interval 



Quantities (|8j)- (|T0|) comprise two type terms differing 
in scaling with respect to r; regular components propor- 
tional to r and anomalous one scaling as s/r. In the 
present section only the latter terms are under consider- 
ation. In deriving the Fokker-Planck equations the divi- 
sion of them by r gives rise to the singularity (t) -1 / 2 . Ex- 
actly their cofactors quantify the influence of the bound- 
ary on Markovian processes and setting them equal to 
zero we can relate the boundary values of the Green 
function G(r,i|ro,io) to the physical properties of the 
medium boundaries. 

Assuming the time scale r to be sufficiently small the 
medium in a certain neighborhood Q s of a boundary 
point s € T is treated as a homogeneous continuum with 
time independent characteristics and the corresponding 
fragment of the boundary T is approximated by a hyper- 
plane. In this case it is naturel to chose the coordinate 
system related to a basis e = ex©w, which, in particular, 
reduces the number of the Green function arguments, 

G(r,x^\r) :=G(r,to+T\{OT,xF},to). 

The system origin was located at the hyperplane T such 
that the vector ro = {Or,x(f} can have only one com- 
ponent xff determining the distance between the point 
Tq and the hyperplane T. Then using the general def- 
initions (jgj)- (fT0"l) of the quantities 9l(r,t,r), il*(r,t,r), 
and £ lJ (r, t, r) the anomalous properties of random walks 
near the boundary T are quantified by their singular com- 
ponents *$1 1 {t,x m ) and *£ ij (r,x M ) scaling as y/r. The 
symbol * is not applied to 9^(r, x M ) because it possesses 
no regular component at all. In other words the desired 
quantities are determined by the following means 

f drG(£,x M \r) = 1 - V\{t,x m ) , (69) 

f drSx i G{r,x M \T) = *il i (T,x M ) + 0(T), (70) 

i / dr5x l 5x 3 G(r,x M \r) = *&(t,x m ) + O(r) , (71) 



~M 



~X M 



„M 



where Sx a = x a and Sx 

In order to calculate these boundary singularities we, 
first, fix the value r and introduce a new time scale 
r a <C r. Then the lattice constructed is Sec. [TV] and 
random walks on it are applied to calculate the desired 
quantities. The advantage of using these lattice random 
walks is due to two reasons. First, the choice of the ba- 
sis b — bx © £>m enables us to simulate the continuous 
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Markovian process as independent random walks along 
the directions parallel to the hyperplane T and along the 
vector bjvf . Second, it becomes possible to ascribe special 
features to the nodes of the boundary layer and in this 
way to simulate some physical properties of the medium 
boundary. In particular, it either can absorb a random 
walker or causes it to migrate extremely fast along the 
boundary within a thin layer. Finally, to restor the con- 
tinuous description the limit r a — > is used. 

The implementation of this approach again is based 
on just mathematical manipulations with the probabil- 
ity function for lattice random walks. So only the final 
results are stated here, referring a reader to Appendix [Bl 
for the proof. 

Proposition 2 Let us consider a Markovian system in 
a homogeneous half-space R M+ bounded by a hyperplane 
T and endowed with the basis b — bx © &m described in 
Proposition [7J 



M-l 

£ 

7=1 



b 7 C 7 + b M ( 



M 



The hyperplane T treated as a physical boundary can ab- 
sorb the system as well as force it to migrate fast along 
T. The diffusion tensor D lJ as well as the drift velocity 
v % at the internal points and v T at the boundary T are 
assumed to be determined in the basis b. It should be 
noted that the boundary drift velocity v T is the velocity 
at which the system had moved outside the boundary if it 
would affected by the same forces. 

The continuous motion of the Markovian system is 
imitated by random walks on the lattice constructed in 
Sec. \IV\ with time step r a . Finally the limit r a — * is 
applied. 

Then, first, the boundary absorption and fast transport 
can be characterized by two kinetic coefficients called the 
surface absorption rate a and the surface diffusion length 
It ascribed directly to the boundary T itself, meaning 
these quantities to be independent of the discretization 
time r a . 

Second, random walks near the hyperplane T exhibit 
anomalous properties reflected in the following singular 
means scaling with the time r as ^pr: 



*ilf(r,C M ) 
*Ufc(r,C M ) 



MM a 

n- 1 / 2 
MM u 



£(r,C M ), 
£(t,C m ), 

^-^ 2 M?-/c(t,c m ), 

S„a-K{r, C M ) 



Ms 



D-^l r V, 



aO a (3 



(72) 
(73) 
(74) 
(75) 



Here the label b notes the basis b used, £ M is the distance 
between the point r and the hyperplane T measured along 
the vector 6m, and the function /C(r, C M ) is specified by 
the integral 



K(r,C M ) 



dz 



■. exp 



(C 



M \2 



4£> m t z 



(76) 



In order to represent these boundary singularities in the 
initial basis e Proposition [1] is applied again. The initial 
basis has been assumed to be of the form e = ex © n with 
the unit normal n to the boundary T directed inwards 
the medium. Let Uy 1 — \\u a p\\ be operator mapping the 
boundary basis ex onto the basis bx. Then transition 
from the coordinates {C"}, Q M of a vector r in the basis 
b to its coordinates {x a },x M in the basis e is specified 
by expressions (|46c[) and (|46d[) using the tensor ii a 3 and 
the diffusion tensor D lJ determined in the initial basis 
e. In the vector form these coordinates are related by 
equality 05). The quantities *il|(r, ( M ) and *£%(t, C m ) 
are obtained by averaging variations of the coordinates 
£*. Thereby, they are a contravariant vector and ten- 
sor, respectively, with the latter being proportional to 
the diffusion tensor written in the basis b and reduced 
to the hyperplane T, namely, the tensor T> a 5 a g. The 
value 5Hfc(r, £ M ) is naturally a scalar. Whence it follows 
directly that 



V\{t,x m )=D^Ho--1C{t,x m ). 



*U 1 (t,x m ) =D 



*£^(t, x m ) = D-^l T T)^ ■ K(r, x M ) . 



-1/2 
MM 

-1/2 



(77) 

JC(t,x m ), (78) 



(79) 



Here the coordinate x and £ are interrelated by for- 
mula (|46d[) and the boundary diffusion tensor J)"' 3 is 
specified by expression (|48|) . Formula (|T8[) can be also 
rewritten in the vector form 



il(r, x M ) = D^l? b + l T v T K(t, x m ) , 



(80) 



where the vector b of boundary singularities is given by 
formula (fT()|) . 



VI. BOUNDARY SINGULARITIES AND THE 
BOUNDARY CONDITIONS 



The obtained expressions (|77)) - ((79^ actually directly 
lead us to the final results. First, they relate the singular 
kinetic coefficients to the diffusion tensor and the physi- 
cal characteristics of the medium boundary. Second, they 
reduce the problem of canceling the singularities inside 
a think layer T T adjacent the boundary T which, nev- 
ertheless, is volumetric before implementing the passage 
to the limit r — > 0. Indeed, since all of terms (|77 p - ([79|) 
depend on the coordinate x M in the normal direction via 
the same function JC(t,x m ) the singularities will be can- 
celed at all the points of the layer T r if it is the case at 
the boundary T. Besides, the structure of the function 
K.(t,x m ), namely, expression (f76|) justifies the adopted 
before assumption that the characteristic thickness of the 
layer T T scales with time as r 1 / 2 . 
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Boundary condition for the backward 
Fokker-Planck equation 



As shown in Sec. Ill Bl the boundary singularities that 
appear in the expansion the Chapman-Kolmogorov equa- 
tion leading to the backward Fokker-Planck equation will 
vanish if equality (fl~8|) holds. At first, in order to lighten 
the perception of results let us consider a rather small 
neighborhood of the point s belonging to the boundary 
T wherein it is actually a hyperplane and chose the basis 
ex © n composed of its hyperplane basis ex(s) and unit 
normal n(s) directed inwards the domain Q. Then sub- 
stituting expressions (|TT|) - (|T5)) into formula (fT5| we im- 
mediately get the conclusion that at the boundary point 
s e T the Green function G(r, t|ro,to) with respect to 
the latter pair of its arguments with ro — > s has to meet 
the condition 

M 

D iM (s, t )Vf G(r, t\s, t ) = a(s, t ) G(r, t|s, t ) 



i=l 



- h{s,t ) 



A/-1 



v%(s,t )V s a G(r,t\s,t ) 



a=l 
M-l 



- J2 S Q/3 (s,to)V^V^G(r,i|s,i ) 

a, 13=1 



(81) 



We note that two last terms in expression (fSTj) describe 
effective motion of the system inside the boundary T 
and has the form of the backward Fokker-Planck opera- 
tor (fl~4"|) with the diffusion tensor D a ^ and drift velocity 
«y whose action is confined to the boundary T. In or- 
der to rewrite this expression for a orthonormal basis of 
general orientation we make use of the definition of the 
boundary singularity vector 6(s, i), expression (|16[) . and 
take into account expression l|48p for the surface tensor 
diffusion. Then introducing the backward Fokker-Planck 
operator acting only within the hyperplane T 



?FP B (s,t ){0} =Zt(Mo) 

■ M M 

W (s, to) V? V* + ]T v\ (s, t ) V? 

-ij=l i=l 



(82) 



where as before the symbol stands for a function acted 
by this operator. Then in the vector invariant form 
the boundary condition for the backward Fokker-Planck 
equation is written as 

6(s, t ) ■ V s G(r, t|s, to) = a(s, t ) G(r, t|s, to) 

-£ F p B (s,to){G(r,i|s,i )} (83) 

which is desired formula. 

In deriving expression (|8 lj) the boundary T was treated 
as a hyperplane, i.e. the Euclidian space of dimension 
(A/ — 1) and its local basis ex was used. To write it again 
in the general form underlining the fact that the operator 



£fp b ac ts in this hyperplane only the tensor notions of 
covariant derivatives are used (see, e.g., Ref. [H|])- In 
these terms the action of the operator ^fp b on the Green 
function taken at the boundary T can be rewritten as 

l F p B {s,t ){G(r,t\s,t )} = ir(Mo) 
M-l 



w?(s > to)G(r,t|s > to) i e 



a=l 
M-l 



- S Q/? (Mo)G(r,i|s,i ) ;Q/3 

<*,P=1 



(84) 



In the given case it is no more that another form of the 
corresponding term in expression (|81[) . However, for a 
nonplanar boundary formula (|84[) holds allowing for the 
boundary curvature, where as expression (|82[) loses the 
curvature effect. Its analysis goes far beyond the scope 
of the present paper, so, here we will just ignore it. 



B. Boundary condition for the forward 
Fokker-Planck equation 

The boundary conditions are obtained in a similar 
way. First, we note that the integrand of expression (|2T)|) 
is similar to the boundary relation (fT5|) within the re- 
placement the test function <f>(r) by the Green function 
G(r, 1 1 ro, to) and action of the operators at the argument 
r in stead of ro. This analogy and the boundary condi- 
tion (fBU]) for the backward Fokker-Planck equation enable 
us to reduce equality ([29]) to the following 

6(s,t) • V s 0(s) = <7(s,t)0(s) -? FPB (s,i){0(s)} (85) 

for an arbitrary boundary point s € T. Since the bound- 
ary part of the backward Fokker-Planck equation acts 
only within the boundary T only the left part of expres- 
sion ([85]) contains the first derivative of the test function 
4>(s) in the direction normal to the boundary T at the 
point s. All the other terms are either the boundary value 
of the function </>(s) itself or its derivatives along the hy- 
perplane T. It justifies the adopted previously statement 
that in the vicinity of T the test function <fi(r) can have 
any boundary value <^(s). 

Then noting that the left-hand side of the condi- 
tion pop is just the combination 

6(s,t)- V s G(s,t|r ,i ) 
the last equability converts expression (|30|) into 



<j) ds(p(s 

T 



M 

E 

i=l 



i/ , '(8)J*{G(B > t|r 0> to)} 1 



= - I ds0(s)cr(s, t) G(s, t\r ,t ) 

T 

+ j> ds£ F p B (s,t){(j)(s)} G(s,i|r ,t ) 

T 



(86) 
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The last term in (|86|) using the divergence integral theo- 
rem for the surfaces is reduced to the form 



ds£ F p B (s,t){cj)(s)} G(s,t\r ,t ) 



dscf)(s) l F p F (s,t){G(s,t\r ,t )} . 



Here the operator ^fp f is the boundary forward Fokker- 
Planck equation 



M 



WM){0}=E V * 



]T V ■ I (s, t)&i (s, t) J - ly(s, t)v\ (s, t ) 
j=i \ ' 

(87) 

where again the symbol stands for the acted function. 
Since the test function <p(s) takes any arbitrary values at 
the boundary T equality ([86]) holds for any point on the 
boundary T, i.e. 

n(s) ■ J {G(s, t\r ,t )} = -a(s, t) G(s, t\r , t ) 

+ £ F P F (s,t){G(s,t\r 07 t )} (88) 

which is the desired boundary condition for the forwards 
Fokker-Planck equation. As it should be the boundary 
condition (|88p can be interpreted in terms of mass con- 
servation; the component of walker flux normal to the 
boundary T is determined by the surface rate of walker 
absorption and the rate of fast surface transport with- 
drawing the walkers from the given boundary point. 



VII. CONCLUSION 

The present paper has developed a technique of deriv- 
ing the boundary conditions for the Fokker-Planck equa- 
tions based on the Chapman-Kolmogorov integral equa- 
tion. The idea of the work is illustrated in Fig. [Jj 

The interest to this problem is partly due to the fol- 
lowing. It is well known that the Fokker-Planck equa- 
tions, forward and backward ones, stem directly from the 
Chapman-Kolmogorov equation under additional two as- 
sumptions, the short time confinement of the correspond- 
ing Markovian process and the local homogeneity of the 
medium. There are rather rigorous techniques of deriving 
them from the integral Chapman-Kolmogorov equation 
based on expanding the latter on short time scales in the 
possible limits. By contrast, the corresponding boundary 
conditions are typically postulated applying to the phys- 
ical meaning of the probability flux and the analogy be- 
tween the forward Fokker-Planck equation and the mass 
conservation law. 



Integral 
Chapman-Kolmogorov 
equation (CKE) 



X 



expansion of CKE in 
two limits for intermediate time 

t, -> l + or S, -» t - 



Forward and backward 
Fokker-Planck 
equations (FPE) 



Standard approach: 

Analogy between 
mass conservation law 
and forward FPE 



Boundary conditions for 




forward FPE 




backward FPE 



FIG. 7: Illustration of the main purpose of the present work 
represented by grey directed line. 



However such simple arguments can fail in dealing with 
more complex Markovian processes like sub- or super- 
diffusion, for which the Fokker-Planck equations with 
fractional derivatives form the governing equations. In 
this case it would be appropriate to have a formal tech- 
nique giving rise to the boundary conditions starting 
from the general description. However, up to now con- 
structing such a technique is a challenging problem. It 
was the case also with respect to the normal Markovian 
processes in continua. 

This paper actually has demonstrated how to do this 
dealing with the normal Markovian processes. The key 
point is the fact that the medium boundary breaks 
down the symmetry of random walks near it. As re- 
sult, the coefficients in the corresponding expansion se- 
ries of the Chapman-Kolmogorov equation are endowed 
with anomalous features called the boundary singulari- 
ties. Namely, they scale on short time scales as St^ 1 / 2 . 
Since the probability distribution on macroscopic scales 
cannot contain such singularities the corresponding co- 
factors in the expressions for the boundary singularities 
should be set equal to zero, leading one to the required 
boundary conditions. In this way we have shown that the 
boundary conditions of the Fokker-Planck equations are 
also the direct consequence of the Chapman-Kolmogorov 
equation supplemented with some rather general assump- 
tions about the properties of the medium boundary. As 
it must the boundary conditions obtained in this way 
match mass conservation. 



APPENDIX A: PROOF OF PROPOSITION Q] 

At the first step the initial basis e of the half-space 
R M+ is assumed to comprise a certain basis 



ex = {ei,e 2 , . . .e M -i} 
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of the hyperplane T and its unit normal n directed in- 
ward R M+ , i.e. e = ex © n. Then the results to be 
obtained will be represented in invariant form where ap- 
propriate, enabling us to write the general expressions. 
The diffusion tensor Z? y is regarded to be determined 
beforehand in the initial basis. Besides, as before, to 
simplify reader perception the Greek letters will be used 
to label tensor indices corresponding to the hyperplane 
T. 

Let us consider a new basis b = by ®b M of the same 
structure except for the last vector b M ; it can be not 
normal to the hyperplane T. An one-to-one map between 
the two bases, e <^4> b, determines a linear transformation 
Li of the space K M mapping, in particular, the hyperplane 
T onto itself. This transformation U = ||£/^-| is specified 
by the relationship between the basis vectors 



M-l 



M-l 



E 



£b 

a=l 



M 



(Al) 



Here the tensor u a p represents an operator Uy acting in 
the hyperplane T whereas the tensor w° (in T) and the 
coefficient uj m ^ complement it to the operator 14, 
namely, 



u p — u f3: 



U M — W J 
M _ , ,M 



(A2a) 

U M p = 0, U M M =u M . (A2b) 

According to the rule of tensor transformations (see, e.g., 
Ref. [23|) in the basis b the diffusion matrix has the com- 
ponents 

M-l 

= u a y y D^' +uj a uj f3 D MM 

7,7 ' =1 (A3a) 

M-l v ; 

7 = 1 

M-l . 

.c ■ 1 Y u % DlM + u a D MM ) , (A3b) 



j^aM _ , ,M 



7=1 



D 



MM _ I^My jjMM 



(A3c) 



Correspondingly, an arbitrary vector x l = {x a ,x M } is 
converted as 



M-l 



x a = Y u V 7 + uj ° lxM ' 
7=1 

X M = LJ M X M . 



(A4a) 
(A4b) 



Currently there is no restrictions imposed on the basis b 
(except for its general structure). Now let us choose a 
specific version of the tensor ui a that eliminates the off- 
diagonal elements of the diffusion tensor in the basis b. 
By virtue of (fAlb)) it is 



M-l 



mm L-i -y 



£)MM 



(A5) 



Here the division by D MM is possible because accord- 
ing to definition (|35|) the diagonal elements of diffusion 
tensor are positive, in particular, D MM > except for 
the case where the system motion along the direction 
n is rigorously deterministic. However, setting further 

D MM _^ + q the latter 

case can be also allowed for. The 
substitution of (|A5|) into (|A3a|) yields 



D 



ap _ 



M-l 

E 

'-1 



7,7 



where the object 



ap _ j~,aP 



D ap - 



JjMM 



jjaM jjPM 



(A6) 



(A7) 



is a tensor within the hyperplane T because, up to now, 
the collections of vectors ex and bx are general bases of 
this hyperplane. 

The tensor D a ^ is symmetric and positive definite. 
The latter property stems directly from inequality ([4*0]) 
written for an arbitrary covector l a of the hyperplane T 
with the component 



hi — — 



1 



M-l 



£)MM 



Y DMll -y . 



(A8) 



7=1 



namely, 



M 



M-l 



D^klj = Y ^ aP hlp > 0. 



(A9) 



a,p=l 



Therefore the basis bx of the hyperplane T can be chosen 
to be orthonormal one wherein the tensor T) a " takes the 
diagonal form with the diagonal components being pos- 
itive values, i.e. D Q ' 3 = = D a p = V a S a p [24]. This 
basis bx is made up of the eigenvectors of the operator 
D := 1 1 55 g || whose eigenvalues are {D a }. For example, in 

the initial basis ex the tensor 2)^ is related to the tensor 
l S a " by the expression 



M-l 

»2 = E 

7=1 



where g a/3 := (e a - ep) 



is the matric tensor of the hyperplane T. 

The choice of the given basis by specifies the transfor- 



mation matrix ||u a ^|| which together with expression (|A5|) 
gives us the vector bu and the corresponding component 
T>m of the diffusion tensor. Namely, first, substituting 
(|A5|) into the latter equality of (|A1[) and taking into ac- 
count the former one we write 

^ M-l 

b M w M = n+— mi Y b a u a J ir' M 



a, 7=1 
M-l 



— V e £> 7M 

MM £^ e ~t U 



JjMM 



(A10) 



7=1 



7=1 
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In the invariant form this expression can be rewritten as 



b M = - e i Dl3 ( e J ■ n ) > 



OJ 

where the normalization factor uj 



(All) 



M 

E 

i,j,k,p—l 



D ik D™{ ei -e ){e k -n){e p -n) 



1/2 



(A12) 



is due to the vector bu being of unit length. Since the 
obtained expressions (|A11[) and (|A12|) are of the tensor 
form and are scaler in this meaning they hold within any 
basis, proving formulae ([IT]) and (j4"2"j) . 

Second, according to (|A10|) and (|A1 1|1 the coefficient 
uj m = lo/D mm . Thereby expressions (|A3cp and (|A12[) 
give us the diffusion tensor component V m related to the 
vector in the basis b 



M-l 



D*i ( ei ■ n)( ej ■ n) 



(A13) 



and written in the invariant form. Formula (|44[) is 
proved. In addition, expressions (|A4[) and (|A5[) together 
with the equality ui = ojmD mm immediately lead to for- 
mulae (|46al) and (|46b[) . 



Finally, we need the transformation Uj- 



of 



the hyperplane T that is inverse to the transformation 
Ur = \\u° U|; its components obey the equality 



M-l 

E 

7=1 



(A14) 



It exists due to the transformation Ut being one-to-one 
map of the bases e and b. Then the inversion of equali- 
ties (|A4[) with u) a given by expression (| A5|) yields formu- 
lae (|46c[) and (|46d[) . Inverting now relationship (|A6[) and 
taking into account the tensor D a P to have the diagonal 
form T> a 5 a /3 in the orthonormal basis b we directly get 

M-l 

7=1 

which together with expressions (|A7|) gives rise to for- 
mula (|47p . The Proposition is proved. 



APPENDIX B: PROOF OF PROPOSITION [2] 

The homogeneous half-space R M+ bounded by the hy- 
perplane T is under consideration and a lattice described 
in Sec. II VI is constructed. It is made up of the node lay- 
ers {Ti} parallel to the hyperplane T with the interplane 
spacing vector a^bM- The individual node arrangement 



of the layers is determined by the vectors of the hy- 
perplane basis bx with spacings {a Q }. In other words, 
the nodes of this lattice are the points 



i\, = 



M-l 

E 

a=l 



n a (a a b a ) + na M bM ■ 



where n is the collection of numbers {nx, n} = {{n a } , n} 
taking any integer value, n a \ l 1 = 0, ±1, ±2, . . ., except 
for the last one; it takes only nonncgative values n = 
0,1,2,... In particular, the points {i" n }T with n = 
form the boundary layer Tq. 



The Markovian process in the half-space 



PA/4 



is sim- 
ulated by random walks on this lattice with hop proba- 
bilities given in Sec. [IV] To find the desired boundary 
singularities we will analyze evolution of the walker dis- 
tribution over the given lattice, i.e. the dynamics of the 
probability 7\ n to find the walker at node n after hop 
t. Here t is the time measured in jump numbers, i.e. in 
units of the hop duration r a . At the initial time t = 
the walker is assumed to be located at a certain internal 
node no. Without lost of generality all the components 
of the index no can be set equal to zero except for the 
last one, i.e. n = {0, 0, . . . , 0, n }. 



Moments of the walker distribution and the 
generation function 



Actually the main purpose of the present appendix is 
to find the zero-th, first, and second order moments of the 
distribution function 7\ m . The zero-th moment quanti- 
fies the trapping effect, whereas the fist and second ones 
characterize the walker propagation in space. Namely, 
the following quantities 



£R„(t,no) =l-EE P *dnr,n}, (Bl) 

n— n-f 

oo 

K(t, n ) = Y, E (" 4 - ™o) Pt,{n Tl n} . (B2) 

n— nx 
^ oo 

^ (t, no) = - E ( ni - <) ( nJ - n o) V - 



n— nx 



t.{nr .n} 

(B3) 



have to be calculated. Here the index i is used as a 
general symbol for one of the indices {a}, M. In order to 
do this the generation function and its analogy written 
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for the boundary nodes only 

oo 

t— ny 

n=0 

(B4) 

OC 

g (s, k T ) = E E e- st+i(kT nT } n, {nT ,o} 

t= rir 

= lim [e- pno G(s,p,k T )l (B5) 

p — >OC 

are introduced, where the complex arguments s, p have 
the positive real parts, Re s, Rep > 0. It should be noted 
that the traps are not included into these sums. The dis- 
crete Laplace transforms of the desired functions (|B1() 
(IB3|) are directly related to the generation function. In- 
deed 



m a (s, n ) = V e~ st m a (t, n ) = 1 - G(s, 0, 0) , 
U {l - e ) 

(B6) 

OO 

K(s,n ) = J2e~ St K(t,n ) = V, G( S ,p,k T )| PikT=0 , 

(B7) 

OO 



n=0 



= 5 ViV 3 -G(a,j»,k T )| 1>ikT=0) (B8) 

where the operator Vj is V Q = —id^a if the index i = a 
is one of the indices of the hyperplane T and V m = —d p 
for the index i = M. 



Master equation for lattice random walks and its 
general solution 

To hnd the generation function for the discrete ran- 
dom walks under consideration the corresponding mas- 
ter equation is applied. For an internal node n = {nx, n} 
with n > 2 it takes the form 



(B9) 



Here prime at the sum denotes the index m running over 
all the nearest neighbors of the given node n and accord- 
ing to expression (|57|) the corresponding hop probabilities 
can be represented as 



1 + CiXi 

2M 



(BIO) 



where e, = T a v l M/cii are some small quantities scaling 

1/2 

with r a as £j tx r a and the value Xi — il stands for 
hops along the basis vector b; or in the opposite direction, 



i.e. the hop to the node with m l = ri 1 ±1 and m? = n J for 
j 7^ i. For the nodes of the layer Ti the master equation 
becomes 



t+l,n 



,Prr 



(Bll) 



Here again prime at the sum has the same meaning ex- 
cept for only internal neighboring nodes being taken into 
account, {rih,n} is the pair of nodes belonging to the 
boundary layer To and the adjacent internal layer Ti 
that are related to each other via walker hops, and the 
hop probability Pi is determined by expression (|60p . The 
walker distribution function P nb ,t in the boundary layer 
obeys the equation 

V t+1 , nb = E Vt, mb PrP^ bnb +V t , n P nnb . (B12) 

We remind that the jumps inside the boundary layer can 
be complex and comprise individually g elementary hops. 

In this case the multihop probability Pm b n b is determined 
by formula (frj5|) . The one- hope probability along the ba- 
sis vector bo, provided the walker remains in the bound- 
ary layer Tq is 



p(l) _ 1 + e a Xa 

mbn " 2(M - 1) ' 



(B13) 



where = T a VyM / a a is again a small parameter scaling 

as oc Tq^ 2 . The values e Q quantify the asymmetry 
of hops in the boundary layer Tq. In particular, these 
complex jumps are characterized by the means 



V n a P {9} - - e T 

4i °" 6 ~ (M - 1) a 



(n a ) r = 
<n-W> T = E « Q ^ P oni 

9 



(B14) 



5 a 



(M-l) " u/3 ^ (Af-1) 2 CqC/3 
Finally, the master equation for the traps is 

VlXU^^+Vt^Ptr. (B16) 

The hop probabilities Pi, Pt r , are given by expressions 
(frjQ)) and the kinetic coefficients of walker jumps inside 
the boundary layer To are specified by expressions (|61[) . 
(frj2")> . and (|65[) . At the initial time the walker distribution 
meets the condition 

7>i=0,n = ^n„ • (B17) 

To solve this system of equations we substitute (IB9[) . 
(|BTTj) . and (|BT2|) into definition JB4]) of the generation 
function G(s,p, kx) and after succeeding mathematical 
manipulations get the following equation (see comments 
about its derivation just after formula (|B20p ) 

[e s -$(p,kT)]G( S ,p,k T ) 

= e s _ e Pno [ $ fa kx ) _ ^ ( p) kx) ] ff ( S; kx) (B1 8) 



5(3-1) x 



(B15) 
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relating the given generation functions G(s,p, kx) and 
<7(s,kx) to each other. Here the following functions 



$(p, k-, 



<f>ip, ki 



M 



(coshp — cm sinhp) 



M-l 



X! ( cos ^* +ie Q sinfc Q ) 

a=l 
(1 ~ <Ta) 



(B19) 



M 

(M - 1) 
M 



-v 



M-l 



^exp[i(k T .n T )]P^ T (B20) 



have been constructed in deriving equation (|B18[) . 



Comments on deriving equation (IB18|) The key fragments 
of deriving equation (|B18[) are outlined below. The con- 
version in (IB4D from t — >• t + 1 leads to the line 



where 



G(s,p,k T ) =e s ^(s,p,k T ) + l, 



0(*,P,k T ) = EE e " S *" P( "" Tl0)+i(kT ' nT) ^+^{nx,n} 



£= riT 

n=0 



and the initial condition (|B17p has been taken into ac- 
count. Equations (|B9)) . (|B11|) . and (|B12|) relating two 
succeeding steps of random walks are substituted into the 
latter expression. As a result the terms in sums (|B9[) 
(|B17p matching the interlayer hops split it into two parts 

g(s,p,k r ) 3>i(p)G(s,p,k r ) 

+ e pno [^i(p)-*i(p)]ff(*,kr) 

with the latter summand caused by that the boundary 
nodes differ from the internal ones in properties. In their 
turn the components of sums (|B9j) - (|B17j) describing tran- 
sitions between a given node n and the nodes of the same 
layer also split the term Q(s,p, kx) into two parts 

g(s,p,k r ) <S> 2 {k T )G{s,p,k T ) 

+ e p "° [0 2 (k T )-$ 2 (k T )] 5 (s,k T ), 

where the latter summand is due to fast diffusion in the 
boundary layer. The combination of the two last lines 
gives equation (|B18p with $ (p, k T ) = $i(p) + $ 2 (kr) 
and 4>{p,k T ) = fa (p) + fa (k T ). □ 

The generation function G(s,p, kx) has no singulari- 
ties in the region Res, Rep > 0. Thereby the left hand- 
side of (|B18|) is equal to zero when e s — $ (p, kx) = 0. 
Resolving the latter equality with respect to the variable 
p we obtain a function p — w{s 1 kx) defined by the equa- 
tion 



$ [oj(s, kx),k T ] 



(B21) 



which specifies the locus in the space {s,p, kx} where 
also the right hand-side of equation (|B18|) has to be equal 
to zero. The latter enables us to write immediately the 
boundary generation function in the form 



ff(s,k-, 



exp [—zu(s, kx)^o] 



1 



'[ro(s,kx),kx] 



(B22) 



Expressions (|B18[) and (|B22|) actually solve the problem 
giving us the following expression for the generation func- 
tion 



G(s,p,k T ) = 



1 



1 - e- s 
1 j>(p,k T )-l] 



[e s -$(p,k T )] [ [l~e- s ] 
p -[»( s ,kr)- P ]n [4> (P, kx) ~ g {p, kx)] 



[l-e-^(w(s,kx) ) k T )] 



(B23) 



where the first summand is the image of the delta func- 
tion 7\ n = <5nn n °t contributing into one of the quanti- 
ties (|Blj) - (|B3|) . the second term is due to random walks 
over the internal nodes, and the last one is caused by the 
boundary effects. Formula (|B23[) specifies the desired 
generation function in the general form. 



Limit of multiple-step random walks on small 
time scales 



In order to find the Laplace transforms (|B6|) (|B8|) it 
suffices to expand the generation function G(s,p, ky) 
into the Taylor series with respect to the arguments p 
and kx with cutting off the series at the second order 
terms. However, in the case under consideration there are 
additional assumptions simplifying essentially obtaining 
the desired results. First, only random walks with many 
steps are of interest because the hop duration r Q has been 
chosen to be much less then the observation time interval 
r of the analyzed Markovian process, r a <C t. It means 
the inequality s C 1 to hold. Second, the time interval r 
is regarded as any small value. So only the components 
of moments (|B1|) — (|B3[) that are characterized by scaling 
T d with the exponent d not exceeding unity, d < 1, are 
to be taken into account. With respect to the generation 
function G(s,p, kx) the latter assumption is converted 
to the statement that all the components of itself and its 
derivatives calculated at the point {kx = 0, p — 0} that 
scale with the argument s as s~ d and have the exponent 
d exceeding two, d > 2, can be ignored. 

At the point {kx =0, p = 0} according to their def- 
inition ljB19l> . (jB2U|) the function $(0,0) = 1 and the 
function 



0(0,0) 



1 

M 



where the coefficient a a is considered to be a small pa- 
rameter, which is justified in the limit r a — > as will be 
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seen below. Thereby in the adopted assumptions expres- 
sion (|B23p for the generation function can be rewritten 



& (s,P, k T ) = - H 5 

+ e - ro(s ,o )na ^(p,k T ) - $(p,k T )] 

s[s + l-0[ro(s,O),O]] v y 

The expansion of the functions $ (p, kx), </> (p, kx) with 
respect to p and kx at the required order is 



$(p,k, 



CMP 

M 



2M 

M-l 



M 2 ( ie< * fea 2 K 



(B25) 



and 



</>(p, k T 

M-l 



Af 



a=l 



1 _ ^ _ JL + J_ 

M Af 2M 

M-l , 
a, /5=l ^ 



g 

2Af 



(B26) 



In deriving expression (|B26p formulae (|B14p . (|B15|) have 
been used. The substitution of the generation function 
written in form (|B24|) with approximations (|B25[) . (|B26|) 
into relations (|B6] ) -([B8 )) yields 



9\ a (s,n ) = -j^fca (s,n ) , 



1 



& a (s,no) = —lC a (s,n 



il«(s,n ) 
Zf (s,n ) 



M 

(.9 - 

Af 

(.9-1) 



CM 

Ms 2 



/C a (s,no) + 



2Af 



£>a£ 



Af - 1 



Ms 2 ' 
/C a (s, no) 



(B27) 
(B28) 
(B29) 



Sat 



2Ms 2 



-i MM 



(s,n ) 



2Ms 2 



(B30) 
(B31) 



the mean £" M (s,no) is equal to zero. Here the function 
fca {s, no) is defined by the expression 



fC a (s,n ) 



exp [—w{s, 0)no] 
s[s + l- (j>(w(s,0),0) 



(B32) 



and we have ignored some insignificant terms where ap- 
propriate. 

Previously in the given appendix we measured time t 
in units of the hop duration r a and spatial coordinates 
{C 1 } in units of the lattice spacings {a^} within the frame 
b. Now let us return to the initial units and deal with the 
corresponding spatial correlations. To do this, first, func- 
tions (|B28[) - (|B31|) should be multiplied by the spacings 



clm and a a , or their products a a ap and a 2 M , respectively. 
Second, the dimensionless Laplace argument s has to be 
replaced by the product sr a , because previously when 
applying to the discrete Laplace transformation the re- 
placement 



st 



t 



has be used obliquely. Third, for further converting 
the discrete Laplace transformation into continuous one 
within the replacement 



t/r a =0 



dt{.. 



all the functions (|B27[) — (|B3ip must be multiplied by the 
time scale r a . 

Leaping ahead we note that the absorption coefficient 
a a has to scale with r a as a a oc */t£. As before noted the 
coefficients {e,;} also behave in this way. Therefore the 
observation time interval r can be chosen to be so small 
that the solution of equation (|B2ip become 



w(sT a ,0) = \J 2Msr a 



(B33) 



and function (|B32p matches a continuous Laplace trans- 
form 



T a K a (sT a ,no) = \j 7r-IC(s, Co) 
2r a 




(B34) 



given by the expression 



K(s,c<r) 



-3/2 



exp 



Co 



s 



(B35) 



with Co = aA/^o being the distance from the node of the 
walker initial position to the medium boundary T along 
the vector &m- 

Indeed, first, if we ignore the second term on the 
right-hand side of expansion (|B25p the solution of equa- 
tion (|B2"T|) for sr a < 1 and k T = is of form (|B33|) . It is 
justified when w 3> cm, which is equivalent to the con- 
dition s 3> v m /T>m or t < T>m/v m . Second, according 
to expansion (|B26|) the denominator in expression (|B32p 
at the leading order is 



[sT a + l-<j>{m(sT a ,0),0)] = 



ZU{STq,0) 

M 



provided zu(sT a ,0) 3> a a - Because a a ~ \j£i~ a , where 
e is some constant, the latter inequality is reduced to 
the following s»£ and r -C e. Since the time interval 
is an arbitrary small value the two inequalities can be 
adopted beforehand. Whence formulae (|B33|) and (|B34|) 
follows immediately for the spacing o,m given by expres- 
sion (|56bp . 
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3. Continuum limit and a 5-boundary model 

To get the final results we analyze the obtained expres- 
sion in the limit r a — > 0. The probability distribution 
T^t.m of the lattice random walks can be treated as the 
discrete implementation of the Green function G(r, ro, t) 
giving the probability density to find a walker at the point 
r at time t provided it was initially at the point ro . Using 
the Green function G(r, r ,t) the means under consider- 
ation are written as the following moments 



5K(t,Co) = l- / dvG(r,r ,t) 

i4(Uo) = 



RM + 

dr(C - Q)G(r,r ,t) , 



(B36) 
(B37) 



dr(C-Co)(C J -C^')G(r,r ,t), 



(B38) 



and their Laplace transforms can be obtained from the 
quantities (|B27| - (|B3l|) in the manner described in the 
previous subsection. As the result we have 



D MM l rVy 1C(S, Co) + -J 



D MM l rT> a K.{s, Co) + 
c2 



(B39) 
(B40) 
(B41) 
S af 3, (B42) 
(B43) 



the component £^ M (s,Co) is equal to zero. Here the 
following characteristics of the medium boundary treated 
as an infinitely thin layer T 



Dmm 
2Mr a 



lr ■■= g 



MDMMTa 



(B44) 



have been introduced and expression (|44|) have been used. 
It should be noted that according to (jB44| the number 
g of elementary hops forming the long distant jumps of 
wallers in the boundary layer To has to grow with r a 
as r a 1//2 in order to retain the effect of boundary fast 
transport in the limit r a — > 0. As a result, the second 
term in the square brackets of expression (|B40j) scales 
as y/r^ because, in turn, the coefficients {ej} vary with 
T a as y/f^. Therefor it vanishes in the limit r a — * 
and the symmetry of the second moments caused by the 
boundary fast diffusion is restored. 

The equality (see, e.g., Ref. [26j |) 



dt 



7rf 



exp 



r 2 



st 



exp 



-Co 



and the Laplace transform of integrals enable us rep- 
resent the inverse Laplace transform K,(t, Co) of func- 
tion (|B35|) in the integral form 



K(t,Co) 



dz 



: exp 



Co 2 1 



AV M t z 



(B45) 



Expression (|B45[1 together with formulae (|B39|) - (|B43|) 
proves Proposition [2] 
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